Codecov Report

Merging #130 into master will decrease coverage by 0.24%.
The diff coverage is 100%.

@@            Coverage Diff             @@
##           master     #130      +/-   ##
==========================================
- Coverage   86.54%   86.29%   -0.25%     
==========================================
  Files          18       18              
  Lines        1390     1394       +4     
==========================================
  Hits         1203     1203              
- Misses        187      191       +4

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update ac7bbee...ca67d46. Read the comment docs.

Hmm, I missed the follow up of @shivin9 in #128, but still I think it is better to not use the KrylovSubspace functions and shoot just for A_mul_B!. CG does Krylov subspace things implicitly, and it feels hacky to initialize a KrylovSubspace of order 1.

I'll do the same bench against #128 soon.

Also: removes the check whether all elements of b are zero. I would say CG should not be responsible for this.

Since there's a test for this, it's fine to change as long as it passes the test. Without that check it used to return all-NaN, but presumably your fix of the starting-point-at-solution bug fixed this.

I think the NaN result is a bug in the current implementation because the initial residual is not checked. With b = zeros(n) and x0 = zeros(n) you start out with a solution, so the iterations should not even start. The current implementation only inspects the residual after the first iteration; this PR does it before the first iteration.

If b = zeros(n) while x0 = rand(n), then CG works just fine.

Removing KrylovSubspace should probably be in a new issue and left intact in this PR. Some solvers in this package already accept KrylovSubspace as input, the intention with making the cg implementation dependant was to enable it as input in the future. So if we are going to remove KrylovSubspace, it should be generalized in the whole package.

If I remember correctly there was an idea to remove KrylovSubspace from the package, right @jiahao? Something along the lines of using it only as a store of vectors and providing normalization functions. I could open a new issue about this. This could also improve readability for people who are more accustomed to see the matrix operations inside the solvers.

@haampie can you leave KrylovSubspace dependency as it was and rebase? To not delay this PR's merge unnecessarily and then we can discuss KrylovSubspace in general for the package.

This PR doesn't remove the KrylovSubspace type, it removes its usage in cg.jl, so this PR can be merged I guess.

Yeah the issue is about removing it's usage in cg, not about removing the type from the package. As I explained before:

Some solvers in this package already accept KrylovSubspace as input, the intention with making the cg implementation dependant was to enable it as input in the future. So if we are going to remove KrylovSubspace, it should be generalized in the whole package.

So the intention was to make cg accept krylovsubspace as input in the future. If we remove it then it can't be done and it should be in a separate issue.

Alright, fixed that

Just for reference, I did a comparison with Trilinos.jl, resulting in the following for CG:

• Your code: 2.7 s
• Trilinos: 1.633 s

I attempted to reproduce your benchmark as accurately as possible here:
https://github.com/barche/Trilinos.jl/blob/master/benchmarks/belos.jl

Thanks @barche for taking the time to do that. It'll serve as a good performance goal to keep in mind! That shows we're already pretty close to Trillinos which is awesome!

Is Trilinos multithreading the sparse matrix vector multiplication? If it is, one could try using something like https://github.com/JuliaSparse/MKLSparse.jl or https://github.com/dcjones/RecursiveSparseBlocks.jl to see if it is the matvec that makes the difference.

@barche did you run Julia with multithreading enabled? In my experience the Julia version runs somewhat faster if you do BLAS.set_num_threads(1) because of all the BLAS level 1 operations. On the other hand @KristofferC might have a good point as well.

@KristofferC No, Trilinos was not multithreading. I'll do some more tests using that when I get back to one of my desktop machines.

@haampie I kept the Julia settings at their default, but switching off threading made no noticeable difference on this 2012 Macbook Air.